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Abstract. We use three-dimensional simulations to study the statistics of 
\ supersonic turbulence in molecular clouds. Our numerical experiments describe 

driven turbulent flows with an isothermal equation of state, Mach numbers 
around 10, and various degrees o f magnetization. We first support the so-called 
1/3- rule of iKritsuk et al.l (|2007af) with our new data from a larger 2048'^ simu- 
lation. We then attempt to extend the 1/3-rule to supersonic MHD turbulence 
and get encouraging preliminary results based on a set of 512'^ simulations. Our 
. , results suggest an interesting new approach to tackle universal scaling relations 

' and intermittency in supersonic MHD turbulence. 

a 
^' 

^ ' 1. Introduction 

H 

' Supersonic turbulence plays an important role in shaping hierarchical internal 

, substructure of molecular clouds (MCs). Since the process of star formation 
begins with the formation of dense cores, turbulence can be responsible for cre- 

. ating initial conditions for star formation. It is hard to access the supersonic 

' regimes typical for MC turbulence in the laboratory and the information avail- 

, able from observations is also limited. Thus numerical simulations currently 

^ I represent the only way to explore the statistics of supersonic turb ulence, the key 

rn 

^ ■ ihle and magnetized interstellar turbulence remained poorly understood. We use 

\ large-scale numerical simulations to shed light on the energy transfer between 

■ scales and on the key spatial correlations of relevant fields in these flows. 

> ' 

^ ' 2. The 1/3-rule for Supersonic Hydrodynamic Turbulence 

H : 

^ In lKritsuk et al.l (|2007al . hereafter K07) we showed that homogeneous isotropic 



ingredient of any successful st atistical theory of star formation (jPadoan et al 



20071 : iMcKee fc Ostrikerl [20071 ). For a long time, the nature of highly 



compress- 



turbulence in a strongly compressible regime typical for the interstellar gas can 
be approximated by modified Kolmogorov's laws since nonlinear advection still 
remains the major physical process at work in the absence of magnetic fields. 
The essence of the required modification boils down to the use of proper density 
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Figure 1. A snapshot of the projected gas density from the 2048'^ Mach 
6 turbulence simulation with PPM. White-blue-yellow colors correspond to 
low-intermediate-high projected density values. 



weighting for the gas velocity, which we call the "1/3-rule" of supersonic turbu- 
lence. Replaci ng the veloci t y u b y p^^^u, one can show that for this quantity 
the 4/5-law of Kolmogorov ( 194ll ) approximately holds at sonic Mach number 
as large as Mg = 6, when the plain velocity power spectrum already scales 
as in Burgers turbulence, i.e. ~ k'^. The first von Karman-Howarth (1938) 
relation between the second-order transverse and longitudinal velocity struc- 
ture functions also still holds approximately at Mg = 6. Although there is no 
rigorous proof for any of these two fundamental exact incompressible laws in su- 
personic regimes, the evidence provided by numerical experiments helps extend 
Kolmogorov's phenomenology to a wider class of problems of interest in turbu- 
lence research and astrophysics. In K07 we showed that the low-order statistics 
of p^/^u are invariant with respect to changes in the Mach number. For instance, 
the slope of the power spectrum is —1.69, and the exponent of the third-order 
structure function S^li) is unity, S-i{£) = {\6u{£)\^) ~ £. 
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Figure 2. Statistics of supersonic turbulence from the 2048^ simulation at 
Ms = 6: (a) power spectrum of projected density (an ensemble average over 
1275 projections); (b) average gas mass M{£) for boxes centered around the 
highest density peaks as a function of box size, the horizontal line corresponds 
to a fractal mass dimension Dm — 2; (c) density pdfs for the 2048'^ simulation 
with a solenoidal driving force (x = 1) and for 1024^ simulation with a hybrid 
driving (x = 0.6, K07); (d) compensated power spectrum of p^/'^u. 



Here we report on our more recent simulation with a grid resolution of 2048^ 
carried out with the ENZO code on Bigben (PSC) and on Ranger (TACC) using 
4096 cores (Fig. [1]). The simulation details can be found in K07, except that 
this new run employs a purely solenoidal large-scale driving force. 

Turbulence statistics derived from these new data with unprecedented dy- 
namic range confirm the convergence of our earlier results obtained with lower 
resolution. Figure [5] provides a brief summary of the new statistics. The power 
spectrum of projected gas density has a slope of —2.01 it 0.01 in the inertial 
range, consistent with the 3D density power index of —1.07 it 0.01 reported 
in K07 (Fig. [2^). The slopes of the power spectra of density and log density 
(not shown) are very close to —1 and to —5/3, respectively. However, it is a 
mere coincidence that at Mach 6 the slopes are represented by these "good" 
numbers, since already at Mg = 10 they become shallower. Our new estimate 
for the "fractal" dimension of the density distribution for the inertial range of 
Mach 6 turbulence. Dm ~ 2.3, agrees well with the previously obtained value 
of 2.4 (Fig. [2)b, see K07 for definition details). In the dissipation range, where 



individual shock fronts are t he domina, r it str u ctures. Dm, = 2 (cf. Kritsuk et al 



2007bl : ISchmidt et aD l2008l : iPan et al.l l2009l : ISchmidt et aD l2009|). Figure [2]c 



compares the probability density functions (pdfs) of the gas density from sim- 
ulations with different driving forces. Both cases agree nicely with lognormal 
distributions at high densities and show some divergence at low densities that 
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can be at least partly attributed to poor statistical representation of violently 
evolving rarefactions due to a limited number of flow snaps hots used to derive 
the pd fs (100 for the 20 48^ data and 170 for K07 data), cf. iLemaster k Stond 
(|2008l 1: iFederrath et~aD (|2008l ^. Finally, Fig. [2]d shows a sample power spec- 
trum of that also confirms the scaling obtained in K07 based on a lower 
resolution simulations. 



3. New Scaling Laws for Supersonic MHD Turbulence 

We carried out a set of pilot simulations of driven isothermal supersonic tur- 
bulence at the sonic Mach number of 10 on 512^ meshes to demonstrate the 



performance of our new PP ML solver for ideal MHD (jPopov Sz Ustvugovl 120071 . 



20081 : lUstvugov et al.l[2009l l. The models are initiated with a uniform magnetic 



field aligned with the x-coordinate direction and cover the transition to turbu- 
lence and its evolution for up to 10 dynamical times. A large-scale {k < 2) 
nonhelical solenoidal force is used to stir the gas in the periodic domain which 
keeps Ms close to 10 during the simulation, see Fig. [3^. The magnetic field 
is not forced and can receive energy only through interaction with the velocity 
field. This random forcing does not generate a mean field, but still leads to am- 
plificatioi i of the small-scale magnetic e nergy via a process known as small-scale 



dynamo (jSchekochihin &: CowlevI 120071 . and references therein). It is implicitly 
assumed that the magnetic Prandtl number Pm ~ 1 in our numerical models. 
Each simulation reaches a steady state with saturated magnetic energy and a 
macroscopically statistically isotropic magnetic field, see Fig. [HJd. The properties 
of this saturated state representing fully developed isotropic compressible MHD 
turbulence is the main focus of this section. While some level of physical un- 
derstanding of small-scale dynamo and its saturation exists, the str ucture of the 



satur ated state is poorly known, even in the incompressible case (jYousef et al 
20071 ). 



Our models are parameterized by the ratios of thermal-to-magnetic pressure 
(3 = Pgas/Pmag- The initial values (3q = 20 and 2 translate into saturated levels 
of the magnetic energy at about 10% and 30% of the kinetic energy, respectively 
(Fig. [3)3). In the /3o = 0.2 case we get a statistical steady state with energy 
equipartition. The three simulations fully cover the transition from highly super- 
Alfvenic regime of turbulence at Pq = 20 with the root mean square (rms) 
Alfvenic Mach number Ma ~ 10 through mildly super-Alfvenic (/3o = 2, Ma_ ~ 
3) to trans- Alfvenic (/3o = 0.2, Ma_ ~ 1), see Fig. [3^. The mean normalized cross- 
helicity cJc = 2 (u • B) /((u^) + (B^)) is contained within ±1.5% for /3o = 20 and 
2, while at /3o = 0.2 the cross- helicity oscillates somewhat more actively with 
(Tc G (—0.073,0.015). Still, to a good approximation, the turbulence remains 
nonhelical. In all three cases PPML keeps the divergence of the magnetic field 
at all times to within |V • B| < 10"^'^. We computed time-average statistics 
over at least four dynamical times for the saturated regimes. The results are 
discussed below. 

The dependence of the density pdf on the value of /3o shows a very clear 
trend for substantially weaker rarefactions in the flows with stronger magnetic 
field (smaller plasma /3), see Fig. ^jp. At a grid resolution of 512^, the aver- 
age minimum density is about W'^-^, lO'^-'^, and lO'^ for (Sq = 20, 2, and 
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0.2, respectively. At the same time, the high-density wing of the pdf preserves 
its lognormal shape and appears insensitive to the intensity level of magnetic 
fluctuations within the studied range of (3o . 

The density power spectra also show some trends with the magnetic field 
strength, although they are less pronounced than in the velocity spectra (see 
below). For example, the power spectrum of the logarithm of projected gas den- 
sity scales as k'^-"^^, k"^'^'^, and k~^-^'^ at /3o = 20, 2, and 0.2, respectively (see 
also Sec. 2 above). The density spectrum for (3o = 20 (nearly nonmagnetized 
medium) at Mg = 10 scales as k~^''^ , i.e. it is shallower than in our nonmagne- 
tized simulations at Ms = 6, where the k~^ scaling was recovered. This result 
confirms K07 prediction that in the limit Mg oo the density spectrum is fiat, 
P{p, k) ~ k^, similar to the white noise spectrum. At a fixed sonic Mach number 
{Ms = 10), as the fiow magnetization increases {Ma drops from 10 to 1), the 
density spectra tend to fiatten, see Fig. [3]d. The same tendency is clear in the 
power spectra of the logarithm of the density, with the slopes around —1.3 at 
Ms = 10. This is consistent with a slope of —1.7 we measured at Ms = 6 and 
Po = oo. We also observe a significant reduction in the extent of the scaling 
interval in the log p spectra as Pq decreases from 20 to 0.2. 

We measured a broad range of the velocity power indices from —1.5 through 
—2 depending on the degree of magnetization (Fig. [3^). As expected, the 
highly super-Alfvenic case Pq = 20 is very similar to K07 results for non- 
magnetized flows, with a Burgers-like scaling of the velocity power spectrum, 
P(u, k) ~ A;"^'^^, and a Kolmogorov-like spectrum for the density-weighted ve- 
locity, P{p^/''^u,k) ~ fc"^ '^ (Fig. EF, also see lKowal fc LazarianI |2007| ) . There is 
a clear trend for the velocity power spectra to get shallower at higher degrees 
of magnetiza tion. We get fc~^' ^ ^ at P n = 2 and k^^'^^ at Pq = 0.2, consis- 
tent with the lLemaster Sz Stond ( 20091 ) measurement for their strong-field case, 
^-1.38 _ Q_Q2 g^j^(j — Q g This result suggests that slopes of the ve- 

locity power spectra around —1. 8 inferred from the observations of molecular 
clouds (e.g., Padoan et al. 20061 ) may indicate a su per-Alfvenic na ture of the 
turbulence there, see however a related discussion in \Li et al.l (|2008l ). 

There is a set of exact scaling laws for homogeneous and isotropic in- 
compressible MHD turbulence analogous to t he 4/5-law of Kolmogorov for or- 
dinary turbulence in neutral fiuids (.Chandrasekharj Il95ll : iPolitano Sz Pouqued 
ll998al M. The MHD laws can be expressed in terms of Elsasser fields, z = 

u±B/V4^ (|Elsasserl[l950l ). as = (^6 zj {e)[6 zf {l)]'^) = -^e±£, where 

(5z||(£) = [z(x + e^) — z(x)] • e, c? is the space dimension, e is a unit vector 
with arbitrary direction, ei is the displacement vector, the brackets denote an 
ensemble average, and summation over repeated indices is implied. Equiva- 
lently, the scaling laws can be rewritten in terms of the basic fields (u, B), 
but in MHD there are no separate exact laws for the velocity or the magnetic 
field alone. It is important not to neglect the correlations between the u and 
B fields or z"*" and z~ fields constrained by the invariance properties of the 
equations in turbulent cascade models. In supersonic turbulence, it is equally 
important to properly account for the density- velocity and density-magnetic 
field correlations. In incompressible MHD, the total energy ((n? + Bf)/2) and 
the cross-helicity (u • B) play the role of ideal invariants and, thus, the (to- 
tal) energy transfer rate = (e^ + e~)/2 and the cross-helicity transfer rate 
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Figure 3. Statistics of MHD turbulence at Ms = 10 and Pq = 20, 2, and 0.2 
from PPML simulations at grid resolution of 512^ points: (a) rms sonic and 
Alfvenic Mach numbers vs. time; (b) mean kinetic and magnetic energy vs. 
time; (c) minimum gas density vs. time; (d) time-average compensated den- 
sity power spectra; (e) time-average compensated velocity power spectra; (f) 
time-average compensated power spectrum of p^/^u for Po — 20; (g) time- 
average compensated power spectra of magnetic energy; (h) time-average 
compensated third-order structure functions 5*3 (£) for generalized Elsaasser 
fields Z*, see text for the definition details. 
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e*-" = (e"*" — e )/2. For vanishing magnetic field B, one recovers the Kolmogorov 
4/5-law, ([(5u || (f)]'^) = —4:/5e£ , where e is the mean rate of kinetic energy transfer 
(|Pohtano Pouauetl[T998bl l. 

Numerical simulations generally confirm these incompressible scalings, al- 
though the Reynolds numbers are perhaps still too small to reproduce the asymp- 
totic linear be havior with a desired precision and the results are sensitive to sta - 
tistical errors (iBiskamp fc Miilleill2000l : iPorter et al]l2002l : iBoldvrev et"alll2006l 1 . 
Often, the absolute value of the longitudinal difference is used and still a linear 
scaling is recovered numerically, while there is no rigorous result for the nor- 
malization constant in this case. The third-order transverse velocity structure 
functions also show a linear scaling in simulations of nonmagnetized turbulence 
(jPorter et al.ll2002. . K07) and the difference of scaling exponents for longitudinal 
and transverse structure functions can se rve as a robust measure o f statistical 
uncertainty of the computed exponents ( Kritsuk fc Norman 20041 ) . In MHD 

simulations, the third-order structure functions of the Elsasser fields (^\SzJ{£)\^^ 

were shown to have an approximately linear scaling too ( Biskamp Sz Miilleill200d : 

Im 

omeni &: Mahdizade h 2008), but again there is no reason for this result to uni- 
versally hold in all situations since the correlations between the fields play 
an important role in nonlinear transfer processes. 

Can the 4/3-law for incompressible MHD turbulence be extended to highly 
compressible supersonic regimes in molecular clouds? As we discussed above, 
a proper density weighting of the velocity, p^^^u, preserves the approximately 
linear scaling of the third-order structure functions at high Mach numbers in 
the nonmagnetized case (for a more involved approach to density weighting 
see iPan et all [20091 ). It is straightforward to redefine the Elsasser fields for 
compressible flows using the 1/3-rule, 7=*= = p^/^ {u±'B / y^Wp) , so that they match 
the original flelds in the incompressible limit and reduce to p^^^u in the limit 
of vanishing B. The new fields can have universal scaling properties in 
homogeneous isotropic turbulent fiows with a broad range of sonic and Alfvenic 
Mach numbers. 

We use data from our MHD simulations to find evidence to support or reject 
this conjecture by computing the third-order structure functions S^^ and 3 
defined above, but now based on the new Z"^ fields. Since we are interested in 
the energy transfer through the inertial interval, we compute the sum of S^^ 

and S^^ which determines the energy transfer rate in the incompressible 
limit. To further reduce statistical errors, we had to combine the transverse 
and longitudinal structure functions, but the absolute value operator was not 
applied to the field differences. The results are plotted in Fig. [3h, which shows 
the compensated scaling for S'3(^) = {S^^ + S^^ + •S'la + Sj^^)/4: averaged over 
about two dozen snapshots for each of the three saturated turbulent states with 
different levels of magnetic fiuctuations. Although the scaling range is rather 
short in 512^ numerical data, each of the three cases clearly demonstrates a 
linear behavior 5*3 (^) ~ i, in contrast to the corresponding scaling of the velocity 
or the magnetic energy which strongly depend on /3o, see Fig. [3^ and g. As a 
consistency check, we computed S'3(^) = ( 5*170 +S'^3)/2 for both plain fields and 



for density-weighted 7=*= fields at /?o = 2. We then compared the scaling S; 
for the two cases and found that the plain Elsasser fields result in a steeper slope. 
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Aa = Oz — ctz ^ 0.23. K07 measured a similar slope difference, Aa ~ 0.31, 
for the third-order longitudinal velocity structure functions in nonmagnetized 
turbulence at Mg = 6. 

While larger simulations are definitely needed to confirm the linear scaling, 
we believe that our preliminary evidence suggests an interesting new approach 
to tackle universal scaling relations and intermittency in compressible MHD 
turbulence. 

4. Conclusions 

We further supported the 1/3-rule of K07 for hydrodynamic supersonic turbu- 
lence with new data from a larger simulation at Mg = 6 on a grid of 2048^ points. 
Based on a series of 512^ MHD turbulence simulations at Mg = 10, we explored 
universal trends in scaling properties of various statistics as a function of the 
magnetic field strength. Our data suggest that the 4/3-law of incompressible 
MHD can be extended to supersonic flows. 
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